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Abstract Future satellite missions dedicated to measuring 
time-variable gravity will need to address the concern of tem- 
poral aliasing errors; i.e., errors due to high-frequency mass 
variations. These errors have been shown to be a limiting 
error source for future missions with improved sensors. One 
method of reducing them is to fly multiple satellite pairs, thus 
increasing the sampling frequency of the mission. While one 
could imagine a system architecture consisting of dozens of 
satellite pairs, this paper explores the more economically fea- 
sible option of optimizing the orbits of two pairs of satellites. 
While the search space for this problem is infinite by nature, 
steps have been made to reduce it via proper assumptions 
regarding some parameters and a large number of numerical 
simulations exploring appropriate ranges for other param- 
eters. A search space originally consisting of 15 variables 
is reduced to two variables with the utmost impact on mis- 
sion performance: the repeat period of both pairs of satellites 
(shown to be near-optimal when they are equal to each other), 
as well as the inclination of one of the satellite pairs (the 
other pair is assumed to be in a polar orbit). To arrive at this 
conclusion, we assume circular orbits, repeat groundtracks 
for both pairs of satellites, a 100-km inter-satellite separa- 
tion distance, and a minimum allowable operational satellite 
altitude of 290 km based on a projected 10-year mission 
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lifetime. Given the scientific objectives of determining time- 
variable hydrology, ice mass variations, and ocean bottom 
pressure signals with higher spatial resolution, we find that 
an optimal architecture consists of a polar pair of satellites 
coupled with a pair inclined at 72°, both in 13 -day repeating 
orbits. This architecture provides a 67% reduction in error 
over one pair of satellites, in addition to reducing the longi- 
tudinal striping to such a level that minimal post-processing 
is required, permitting a substantial increase in the spatial 
resolution of the gravity field products. It should be empha- 
sized that given different sets of scientific objectives for the 
mission, or a different minimum allowable satellite altitude, 
different architectures might be selected. 

Keywords Time variable gravity • GRACE • Temporal 
aliasing errors • Constellations • Satellite geodesy 


1 Introduction 

The Gravity Recovery and Climate Experiment (GRACE) 
satellite mission has demonstrated the ability to measure 
mass variations on the Earth at monthly to sub-monthly res- 
olution (Tapley et al. 2004a; Bruinsma et al. 2010). There is 
tremendous potential in the science applications of this data 
set, including quantifying variations in continental hydrol- 
ogy, measuring ocean currents and ocean bottom pressure 
variations, and sensing the loss and accumulation of ice in the 
polar regions as well as glaciated regions around the globe. 
Because GRACE senses these signals at large spatial scales 
(~ 400 km), mass variations in smaller drainage basins, for 
example, cannot be quantified by GRACE. Thus, it is of pri- 
mary interest for future missions to have increased spatial 
resolution in the time-variable gravity field products. One 
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potential benefit of this data set would be to aid in regional 
water management practices. 

In designing future missions, we have now reached a point 
where the spatial resolution will no longer be limited by 
sensor errors, but rather by temporal aliasing errors (Wiese 
et al. 2011; Loomis 2009; Visser et al. 2010), assuming that 
future missions will use a laser interferometer for the inter- 
satellite ranging (Pierce et al. 2008; Bender et al. 2003), and 
the spacecraft will fly drag-free at lower altitudes, similar to 
the Gravity field and steady-state Ocean Circulation Explorer 
(GOCE) (Drinkwater et al. 2007) mission. Temporal aliasing 
errors are defined as errors due to temporal undersampling of 
geophysical signals of interest (e.g., continental hydrology) 
as well as mis-modeling unwanted high-frequency mass vari- 
ations (e.g., atmosphere and ocean tides). 

We aim to directly reduce the effect of temporal aliasing 
errors by adding a second pair of satellites, thus increasing 
the sampling frequency of the mission and improving the 
spatial resolution of the derived gravity field products. In 
theory, if enough pairs of satellites were placed in proper 
orbits, one could sample the gravity field at a high enough 
frequency such that temporal aliasing errors would be largely 
eliminated. However, there could still be a problem of signal 
separation given this scenario. Having a dozen, if not more, 
satellite pairs to accomplish such a feat is cost-prohibitive at 
this point. As such, this paper focuses on the more interesting 
question of optimizing the orbits of two pairs of satellites. 

One reason that the GRACE gravity field products have 
limited spatial resolution is the fact that the polar orbiting 
two-satellite collinear architecture has little East- West sen- 
sitivity to the gravity field. This leads to a strong correlation 
between coefficients of a fixed order and the same parity 
of degree, manifesting as longitudinal stripes in the gravity 
solutions. One common method used to remove the stripes 
is to apply a decorrelation filter (Swenson and Wahr 2006) 
to the data (known as ‘destriping’). The gravity fields are 
also typically smoothed with a Gaussian filter (Wahr and 
Molenaar 1998) to average through any remaining errors. 
The resulting product is a map of global mass variations at 
large spatial scales. While substantially reducing the level 
of errors in the gravity solutions, these two processes have 
the undesired effects of removing real geophysical signal as 
well as reducing the spatial resolution of the solutions, par- 
ticularly in regions with high signal to error ratios. Hence, 
it would be advantageous if future missions did not have to 
rely on these techniques. 

While simulation studies have shown that alternate forma- 
tions, such as cartwheels and pendulums, reduce the level of 
striping in the gravity solutions (Elsaka 2010; Wiese et al. 
2009), these types of formations are much more difficult 
to implement from an engineering standpoint. Alternately, 
Bender et al. (2008) suggested that having two pairs of col- 
linear satellites, one pair in a polar orbit and one in a lower 


inclined orbit, would improve the East-West sensitivity to 
the gravity field variations, thus leading to a reduction in the 
level of striping. This concept is explored in this paper. 


2 Orbit design considerations 

The goal of this paper is to optimize the orbits of two pairs of 
collinear satellites for improving the spatial resolution of the 
derived gravity field products for future missions over what 
GRACE provides. The search space for such a problem is 
extremely large, and is further complicated when considering 
that the selected orbits will be a strong function of the sci- 
ence goals of the mission. For instance, if the primary goal of 
the mission is to determine continental hydrology (excluding 
ice) at small spatial scales, then one might place the satellites 
in orbits with dense coverage over these regions, but less cov- 
erage over the polar regions. This would most likely result 
in decreased sensitivity to determining ice mass variations in 
Greenland and Antarctica. However, if the primary science 
objective is to determine ice mass variations in Greenland, 
then a different mission architecture would be selected. This 
study assumes that the science goals of the mission are to 
determine continental hydrology, ice mass variations, and 
ocean bottom pressure signals over the entire globe, with 
each area of science being weighted equally. 

Considering strictly the satellite orbits, one can character- 
ize the mission performance, P, given two pairs of collinear 
satellites, via the following: 

P = f(X u X 2 , Au 1 ,Au 2 ,L). (1) 

In Eq. 1, X\ and X 2 are the state (position and velocity) 
of the lead spacecraft of the first and second pair of satel- 
lites, respectively, Aiq and Av 2 are the separation distances 
between the first and second pairs of satellites, respectively, 
and L is the amount of time that data are collected to form a 
single solution. It is most convenient to represent the state of 
the spacecraft in terms of mean Keplerian orbital elements 
as 

it, £2i, cut, vi) 

X 2 — f ( a 2 , e 2 , i 2 , ^ 2 , ® 2 , V 2 ) 

Here, a is the semimajor axis, e is the eccentricity, i is 
the inclination, £2 is the longitude of ascending node, co is 
the argument of perigee, and v is the true anomaly. Coupling 
Eqs. 1 and 2, one can see that the mission performance of 
this type of architecture will be directly related to 15 param- 
eters. Adding additional satellite pairs will increase the num- 
ber of variables by seven for each pair of satellites added. It 
is desirable to reduce the number of independent variables 
and narrow down the search space by making appropriate 
assumptions. 
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First, the inter-satellite separation distances, defined as 
A i; i and A V 2 , will likely be chosen based on the satellite-to- 
satellite ranging instrument requirements. Future missions 
are likely to use a laser interferometer, for which a 100-krn 
separation distance is chosen as a trade-off between instru- 
ment performance as well as relative accuracy in determining 
short wavelength and long wavelength features in the gravity 
field (Wiese et al. 2009). Fixing this distance allows us to 
eliminate two of the variables, Aui and Ai’ 2 , from the search 
space. 

Next, it can be assumed that the spacecraft should fly in 
circular orbits to minimize any relative changes in distance 
due to having eccentric orbits, as GRACE does. Fixing the 
eccentricity to zero eliminates two additional parameters. 

Given circular orbits, the argument of perigee, co, becomes 
ill-defined. Hence, we can now define the argument of lati- 
tude, u, which is equal to the sum of the true anomaly and the 
argument of perigee (n = w + v). The argument of latitude 
defines the position of the satellite in its orbit about the Earth 
with respect to Q. While this is a parameter that could have 
an effect on the gravity solution, it is impossible to deter- 
mine what the optimal satellite position should be due to the 
extremely complex nature of the problem. For example, it 
would be optimal if, during a flooding event, the satellite 
flew over the region of interest. However, it is impossible to 
know when this event might occur in the future, making it 
very difficult to optimize. 

While optimizing u\ and ui independently is not feasi- 
ble, one could think of optimizing the relative difference in 
the argument of latitude between the two pairs of satellites 
in an effort to meet certain temporal groundtrack crossing 
constraints (i.e., the second satellite pair will fly over a loca- 
tion on the Earth a specified amount of time after the first 
satellite pair flew over the same location). The same argu- 
ment holds for the longitude of ascending node, Q, in a spa- 
tial sense. That is, in an absolute sense it is impossible to 
determine what the optimal values for and Fb should 
be, since we cannot predict the time and location of mass 
variations on the Earth years in advance. However, the rela- 
tive difference between the ascending nodes of the two pairs 
of spacecraft could be optimized to provide a required spa- 
tial constraint on the combined groundtrack pattern of the 
two satellite pairs. Thus, Q i and along with u \ and uj 
can be reduced to two new parameters: Af 2 i 2 , and Am 12. 
The first, AG 12 , provides a spatial constraint on the ground- 
track pattern of the two satellite pairs while the second, Am 12 
provides a temporal constraint on the groundtrack pattern. 
It is expected that Am 12 can only be optimized if the peri- 
ods of both satellites pairs are equal to each other, which 
would require that a 1 = aj- Otherwise, there will be a sec- 
ular drift rate in the time that the two pairs of satellites 
cross the same location on the Earth which cannot be 
controlled. 


Next, we can consider the inclination of the satellites. In 
order to provide global coverage of the Earth, at least one of 
the pairs of satellites must be in a near-polar orbit. Thus, this 
can be set as a constraint. The inclination of the second pair 
of satellites, however, is free to vary. The problem has now 
been reduced from one with 1 5 parameters to one with only 
six, and can be represented via the following: 

P = /(fli, a2, h, Am 12 , A£2 12 , L). (3) 

Let us now discuss the semimajor axis of the two pairs of 
satellites. GRACE was launched into an altitude of approx- 
imately 490 km (Tapley et al. 2004b) in which the altitude 
of the satellites continuously decays primarily due to atmo- 
spheric drag forces. Currently, GRACE is at a nominal mea- 
surement altitude of 460 km (Center for Space Research 
2011). GOCE, alternatively, which is designed to measure 
the static gravity field of the Earth to high spatial resolu- 
tion, flies at a much lower altitude than GRACE, at approx- 
imately 255 km. It employs a single-axis drag-free system 
where non-conservative forces are measured in real time and 
are compensated for using thrusters to maintain a nominal 
altitude. GOCE was designed for a lifetime of 2 years and 
is limited by the amount of propellant onboard to keep the 
spacecraft at the nominal measurement altitude. This can be 
contrasted with the GRACE mission, which is currently in its 
ninth year of operations. It should be noted, however, that due 
to the extended solar minimum at the end of the past decade, 
GOCE is now expected to continue performing several years 
beyond its mission lifetime (Fehringer et al. 2010). 

It is envisioned that future GRACE-type missions will 
also employ drag-free technology, allowing one to fly at a 
lower altitude with increased sensitivity to short wavelength 
features in the gravity field. With this in mind, one can set 
a minimum bound on the altitude of the spacecraft which 
depends on many factors, including, but not limited to the 
design lifetime of the satellites, the amount of propellant 
available, the type of thrusters used, the cross-sectional area 
of the spacecraft, and the magnitude of the atmospheric den- 
sity. Some work has been done to this end, by Marchetti 
et al. (2008) and St. Rock et al. (2006), examining the per- 
formance of drag-free control systems in low-Earth orbit, 
and the mission lifetimes associated with various thrusters. 
Figure 1 depicts the results from each respective paper, along 
with the initial estimate for the GOCE mission, assuming the 
same initial mass propellant fraction as the GRACE mission 
(0.18). Note that the results from St. Rock et al. (2006) have 
been scaled down by a factor of two to account for variable 
specific impulse and control system use that was not consid- 
ered in the analysis. 

Given that GRACE is quickly approaching its tenth year of 
operations, it seems both reasonable and desirable to design 
future missions with a lifetime of at least 10 years. Fig- 
ure 1 illustrates that a 290-km altitude allows the satellites 
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Fig. 1 Mission lifetime as a function of altitude assuming an initial 
mass propellant fraction of 0. 1 8 


to remain in orbit for 10 years; thus, this was selected as the 
minimum altitude for this study. Note that this calculation is 
very approximate, and a rigorous analysis of a control sys- 
tem in the appropriate environment would need to be made to 
refine the targeted altitude; however, it is valid as a first-order 
approximation and sufficient for the purposes of this study. 

The last parameter which needs to be discussed from Eq. 3 
is L , the length of time that data are collected. The parameter 
L will depend primarily on the targeted spatial and tempo- 
ral resolution of the mission. In principle, the product of 
the spatial resolution and temporal resolution of a mission 
is constant; that is, given a fixed number of satellite pairs, 
one cannot improve the spatial resolution without sacrificing 
temporal resolution and vice-versa (Visser et al. 2010). Given 
homogeneous groundtrack spacing, the spatial resolution of a 
mission can be approximated by the Colombo-Nyquist rule, 
which states the maximum resolvable degree of the gravity 
field is equal to half of the number of orbital revolutions of 
the satellites (Colombo 1984), guaranteeing gravity solutions 
that are homogeneous in longitude (Visser et al. 2011). While 
larger values of L theoretically lead to better spatial resolu- 
tion, they also allow for greater accumulation of temporal 
aliasing errors. Varying the L parameter should lend insight 
into proper trade-offs between increasing the spatial resolu- 
tion of the solutions and mitigating the effect of temporal 
aliasing errors. 

There is one more point to be mentioned concerning L. 
One drawback of GRACE is the lack of an altitude control 
system. This leads to variability in the groundtrack pattern 
of GRACE, and subsequent variability in the quality of the 
monthly solutions. This was discussed in Klokocmk et al. 
(2008) and Wagner et al. (2006), showing the degradation 
in gravity solutions from GRACE in the fall of 2004 when 
the satellites passed through a 61/4 resonance orbit. Due to 
the success of the GOCE mission implementing a drag-free 
system and maintaining a constant orbital altitude, it seems 
advantageous to consider only repeat groundtracks for future 
missions. Imposing this constraint assures consistent quality 


in the time-variable gravity solutions. It should be noted, 
however, that maintaining repeat groundtracks with a pair of 
satellites will be more demanding than for a single satellite, 
like GOCE. Tolerances on the formation dynamics will be 
set based on instrument performance criteria and how tightly 
one needs to maintain the repeat groundtrack (presumably 
this will be a function of the selected orbit and what other 
repeat orbits exist close to the operational altitude). These 
parameters will need to be given careful consideration should 
such a mission be flown. 

From Kaula (1966), even zonal coefficients in addition to 
a nonlinear (7 2 ) 2 contribute to the secular rate of the node, G. 
However, all terms are second order when compared with the 
contribution of J 2 ; thus, we design the repeat groundtracks by 
considering only this term. Given a desired eccentricity and 
inclination, there are only certain values for the semimajor 
axis which will satisfy the conditions for a repeating ground- 
track. One can obtain the appropriate values for semimajor 
axis, a, by solving the following equation (Vallado 2001): 

C 2 a 7/2 + Cia 2 + C 0 = 0, (4) 


where 

/ 

C 2 = 7 &>e 
k 

1 ' 1 — 


(5) 


Co= 4? 

in which 


2 cos ; + 1 — 5 cos 2 i — ^3 


cos 2 / — 1^ e 


a = 3^fjlJ2Rl 

e=(l -e 2 ) 2 - (6) 


In these sets of equations, /i is the gravitational constant 
of the Earth, J 2 is the negative of the unnormalized C 20 coef- 
ficient describing the oblateness of the Earth, R e is the radius 
of the Earth, co e is the rotation rate of the Earth, k is the desired 
number of nodal days it takes for the satellites to repeat, and 
/ is the number of orbital revolutions the satellites perform in 
k nodal days. It should be noted that k/ 1 must be irreducible, 
and / is given by 


k (a> + M) 
(u>e G) 


(7) 


In Eq. 7, co and M are the secular drift rates of the argu- 
ment of perigee and mean anomaly, respectively, due to the 
oblateness of the Earth. 

Thus, selecting a particular value for L inadvertently 
imposes an additional constraint on either a 1 or « 2 : that the 
value for a must put the satellite in a repeat orbit. It is not 
imperative that both satellite pairs have a value of k equal to 
that of L, but one pair must. It has been pointed out by Bender 
et al. (2008) that perhaps the most effective way to design the 
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Table 1 Constraints imposed on design criteria to reduce the search 
space for an optimal architecture 


Parameter 

Constraint 

a\, ao 

290 km minimum; repeat groundtrack 

ei, ei 

0 (circular orbits) 

h 

90° 

Avi , Ai’2 

100 km 

L 

L = k 2 ; k\ < £2 


architecture would be to have a lower inclined pair in a longer 
repeat period (RP) coupled with a polar pair of satellites in 
a shorter RP. This would lead to more homogeneous spac- 
ing in the combined groundtrack pattern of the two pairs of 
satellites, since, by nature, groundtracks are more dense over 
the poles than the equator. In this scenario, the lower inclined 
pair would be selected to have a value of k equal to L, while 
the RP of the polar pair of satellites would be allowed to vary, 
but would be constrained to be less than or equal to that of 
the lower inclined pair. Thus, all such combinations should 
be explored. 

Finally, taking into consideration the above discussion, 
Eq. 3 can be rewritten as 

P = f(ki, k2, h, Ami 2, A£2i2). (8) 

In Eq. 8, k\ is the RP of the polar pair of satellites and kj is 
the RP of the other pair of satellites for which the inclina- 
tion can vary. Note that the additional constraints which are 
imposed are that ki = L and that k\ < kj. One additional 
caveat that should be mentioned is that k is expressed in units 
of nodal days, while L is typically expressed in units of solar 
days, since the data processing is usually set up to handle 
daily batches of data. This means that typically a solution 
will have slightly more data (a few hours) than what is taken 
during the full repeat period of the satellites. Table 1 is a list 
of all constraints that were imposed to arrive at Eq. 8. 

Using the constraints listed in Table 1, Eq. 8 has been 
reduced from one that initially was a function of 15 variables 
to one that is now a function of only five variables. Further- 
more, it is expected that the values selected for k\ , kj, and 
(2 will have the most influence on how well the mission per- 
forms. Ah 12 and AQ|o are expected to have much smaller 
impacts. 

It should further be stressed that this type of analysis is 
considerably biased towards the minimum altitude chosen, 
in this case, 290 km. To illustrate this, Fig. 2 shows the clos- 
est altitude to 290 km (without going below it) for different 
values of k\ for a polar pair of satellites. 

Figure 2 illustrates how results could be biased towards the 
minimum allowable altitude. For example, the closest 8-day 
RP groundtrack to 290 km exists at an altitude of 29 1 km, ver- 
sus 374 km for a 12-day RP. The lower altitude given by the 
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Fig. 3 Necessary altitude to maintain a 17-day repeat period at differ- 
ent inclinations 




8-day RP orbit could trump any benefit that collecting data 
for 12 days versus 8 days might add. 

The same analysis can be done if one fixes a value for 
kj, but lets zb vary. Figure 3 shows how the altitude neces- 
sary to maintain a 17-day RP orbit changes as a function of 
inclination. 

Hypothetically speaking, if a 70° inclination were an opti- 
mal value for ij speaking strictly in terms of inclination. Fig. 3 
suggests that this study might find that a 7 1 ° inclination is the 
optimal value for ij since this orbit is 16 km lower in altitude 
than the orbit with a 70° inclination. This analysis shows how 
the altitude, repeat period, and inclination are inherently cou- 
pled together, and an optimal set of orbital parameters will 
be a strong function of the mission constraints. 


3 Methodology 

Numerical simulations were used to explore the full search 
space of parameters for the investigation. This was deemed 
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necessary due to the variable spatio-temporal sampling char- 
acteristics of different mission architectures along with the 
unpredictable nature of geophysical signals and temporal ali- 
asing errors. The primary parameters of interest in Eq. 8 are 
and z 2 . The following are the ranges of values explored 
for these parameters: 

ki,k 2 = [4,23] days; h < k 2 (9) 

i 2 = [25°, 90°] (10) 

Using a step size of 5° for z' 2 , it is readily seen that 2,940 
simulations were run to explore the entire search space of 
the k \ , A' 2 - and z 2 variables. After the initial results narrowed 
down a more appropriate range of values for z 2 , a 1° step 
size was used. While a more systematic method could have 
been implemented to explore the search space and reduce 
computation time, this would have led to greater uncertainty 
in the validity of our analysis due to the variable spatio- 
temporal characteristics of geophysical signals, temporal 
aliasing errors, and the sampling behavior of the different 
architectures. 

3.1 Simulation procedure 

The GEODYN software package (Pavlis et al. 2010) (an orbit 
determination program) coupled with SOLVE (Ullman 1997) 
(a linear systems solver), both from NASA Goddard Space 
Flight Center, was used to perform the numerical simulations. 
In an effort to create as realistic of a simulation environment 
as possible, we model the dominant mass variations in the 
Earth system, including the atmosphere, oceans (mass redis- 
tribution due to atmospheric forcing), tides (ocean, atmo- 
sphere, and solid Earth), continental hydrology, and ice 
sheets. Table 2 shows the simulation design. Throughout 
the paper, when we refer to recovering hydrology signals, 
we actually mean recovering continental hydrology, while 
recovering ice signals refers to recovering ice mass varia- 
tions in Greenland and Antarctica. 

The static gravity field model used, EIGEN-GL04C, 
includes data from the GRACE and LAGEOS missions along 
with surface gravity data from altimetry over the oceans 
and gravimetry over the continents (Forste et al. 2008). The 


Table 2 Simulation and model definitions 


Models 

Truth 

Nominal 

Static gravity field 

EIGEN-GL04C 

EIGEN-GL04C 

Ocean tide model 

FES2004 

GOT00 

Atmospheric model 

ECMWF 

NCEP 

Ocean model 

OMCT 

MOG-2D 

Hydrological model 

GLDAS 

None 

Ice model 

ESA 

None 


atmospheric models used are 3-h ECMWF surface pressure 
fields (Klinker et al. 2000) and 6-h NCEP Reanalysis fields 
(Kalnay et al. 1996) while the ocean models are the baroclin- 
ic OMCT model which is used as a dealiasing product for 
GRACE (Flechtner 2007), along with the barotropic MOG- 
2D model, both of which are forced by ECMWF surface 
pressure (Carrere and Lyard 2003). The tide models used are 
of two generations. GOTOO.2 uses almost 8 years of TOPEX 
and Poseidon data, supplemented in shallow seas and in polar 
seas (latitudes above 66°) by ERS-1 and ERS-2 data and 
uses FES94 as an apriori model (Ray 1999). FES2004 is a 
hydrodynamic model that not only assimilates more recent 
altimetry data, but also includes that from TOPEX (Lyard 
et al. 2006). The 3-h GLDAS/Noah land-surface model is 
used to describe continental hydrology (Rodell et al. 2004) 
and the 6-h ESA ice model describes the ice sheet dynam- 
ics in Greenland and Antarctica, and was provided by the 
European Space Agency (ESA) (van Dam et al. 2008). All 
simulations are arbitrarily run for January of 2003. The ice 
model is defined from 1995 to 2006 with 1995 being the ref- 
erence time; as such, the magnitude of the ice signal in 2003 is 
significant. All models are represented to spherical harmonic 
degree and order 100 in the simulations with the exception 
of the NCEP and MOG-2D models, which are represented 
to spherical harmonic degree and order 72. Note that solid 
Earth and atmospheric tides are included in the simulation, 
but errors are not considered. 

Table 2 shows that this particular simulation is designed to 
recover hydrology and ice mass variations in the presence of 
model errors from ocean tides and atmospheric and oceanic 
mass variations, defined as the differences between the two 
sets of models. One can also calculate ocean bottom pressure 
variations from the above simulation by treating the NCEP 
and MOG-2D models as forward models and examining dif- 
ferences with respect to the truth models over the oceans. All 
models are interpolated linearly in time during the simulation 
process. 

The spacecraft are assumed to be equipped with a laser 
interferometer as a replacement for the microwave ranging 
instrument onboard the GRACE mission. The noise assumed 
for the laser interferometer is 5 nm/VHz (Young et al. 1999; 
Mueller et al. 2005; Wiese et al. 2009). Additionally, the 
spacecraft are assumed to fly drag-free, such that they main- 
tain their nominal orbital altitudes. The noise on the drag- 
free system is taken to be 0.01 nm/s 2 /VHz, which is 
similar to that of the GOCE mission. Spacecraft position 
errors are added as white noise with a 1 -cm RMS in the radial, 
alongtrack, and crosstrack directions. Simulations show that 
with these levels of error, temporal aliasing errors from mis- 
modeling atmosphere and ocean signals dominate the error 
budget by more than two orders of magnitude over the level 
of error from the drag-free system and laser interferometer 
(Wiese et al. 2011). Thus, any improvements that certain 
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architectures offer over other architectures will be attributed 
to reducing temporal aliasing errors. 

The simulations are carried out in two steps. First, the 
selected satellite orbits are integrated through both the truth 
and nominal sets of force models independently. Noise on the 
laser interferometer, drag-free system, and satellite positions 
are added to the truth set of measurements. The residuals 
between the inter-satellite range-rate measurements 
coupled with the satellite position measurements are then 
used to estimate a correction to the spacecraft state. This step 
is performed to account for a change in energy of the sys- 
tem due to a different set of force models. The second step 
involves integrating the spacecraft orbits through the nom- 
inal set of force models using the updated spacecraft state. 
The range-rate residuals are then used exclusively to estimate 
any additional corrections to the state of the spacecraft along 
with corrections to the geopotential coefficients. Spacecraft 
position measurements are not used in this step, as the space- 
craft states are converted to baseline elements via Rowlands 
et al. (2002), and nine out of the twelve parameters are con- 
strained during the estimation process. Data are collected in 
one-day batches and SOLVE is used to combine ki days of 
data together to arrive at the final multi-day solution. 

The computation time associated with such a large matrix 
of simulations is expensive, and increasing the degree of esti- 
mation exponentially increases the processing time. As such, 
a subset of simulations was carried out to both degree 60 and 
degree 100 to ensure consistency between the two, in hopes 
of being able to run the matrix of simulations to degree 60. 
It was expected that the two would correlate; however, the 
results were surprising, showing smaller correlations than 
expected. This is discussed in detail in Sect. 4.1. As a result, 
all simulations were run to degree and order 100. 

3.2 Performance metrics 

There are several ways one can quantify the performance of 
a gravity recovery satellite mission, given by P in Eq. 8. As 
discussed in Sect. 2, P depends substantially on what the sci- 
entific goals of the mission are. For this study, we took the 
liberty of defining the scientific goals of the mission to be 
increasing the spatial resolution of the recovered hydrology, 
ice, and ocean bottom pressure signals over what GRACE 
provides. Each area of science is weighted equally; hence, 
one hopes to minimize the error, E, given by 

_ £(H) + £(I) + E (O) 

3 C 

InEq. 1 1 , E (H) represents the error in determining hydrol- 
ogy, E (I) is the error in determining ice mass variations, and 
E (O) is the error in determining ocean bottom pressure sig- 
nals. There are many methods and tools which one can use 
to analyze error and quantify E (H), E (I), and E (O) on both 


global and regional scales. It was pointed out in Han and 
Ditmar (2008) that global error metrics lack sophistication 
as they disregard the different spatial distributions of signals 
and errors, making it desirable to analyze errors on a regional 
basis. 

There are several types of localized analysis one can use, 
including, but not limited to, localized averaging kernels 
(Swenson and Wahr 2002) and a spatiospectral localization 
technique (Han and Simons 2008; Han and Ditmar 2008; 
Simons et al. 2006; Wieczorek and Simons 2005). While 
these are effective techniques to analyze regional mass vari- 
ations, a complete analysis of this kind requires one to select 
multiple regions of varying size, signal strength, error charac- 
teristics, etc., to properly analyze the results. Since this study 
is already computationally expensive, it is desirable to use a 
global metric to quantify E to narrow down the search space 
and identify a select few mission architectures for further 
analysis on a regional scale. This type of regional analysis 
will be addressed in future work. 

In a global sense, one of the most common ways to look at 
error is using degree variances. While this is a valuable tool to 
examine errors in the spectral domain globally, the end user 
of the GRACE data is more interested in what is happening 
in the spatial domain. Thus, it was determined that the metric 
used for calculating £(H), EG), and E( O) would be differ- 
encing the truth signals and the recovered signals to obtain a 
spatial plot of errors. The spatial plot of errors is represented 
on a 1° x 1° grid, and an area- weighted RMS of the errors 
for each signal is calculated and substituted into Eq. 1 1 for 
Zt(H), EG), and E( O). While this is not a perfect represen- 
tation for the performance of a mission by itself, it does give 
a very good indication of how changing k \ , ki, and h affects 
the ability of the satellites to recover the geophysical signals 
that we are interested in. 

Figure 4 illustrates how E is calculated, showing the truth 
signals (left), recovered signals (middle), and the recovered 
signals after being destriped via Swenson and Wahr (2006) 
and smoothed with a 300-km Gaussian averaging radius 
(right), for recovering both hydrology and ice mass variations 
(top), as well as ocean bottom pressure signals (bottom). This 
simulation is for a single pair of polar orbiting satellites in 
a 13-day RP at 299 km. The figures have been truncated at 
degree 60 and are expressed in cm of equivalent water height 
(EWH). Note that in this case, the truth signal for hydrology 
and ice is defined as the 13-day average of GLDAS + ESA, 
while the truth ocean signal is defined as the 13 -day aver- 
age of ECMWF + OMCT. Any non-zero mean between the 
atmosphere and ocean models over the time span of interest 
is compensated for in the post-processing when estimating 
hydrology and ice signals. 

When calculating E, one uses a spatial plot of the errors 
which is obtained by differencing the truth signals, given by 
the left set of figures, with the recovered signals, given by the 


■£) Springer 


D. N. Wiese et al. 




Fig. 4 Truth signals (left), recovered signals (middle), and recovered pair of satellites at 299 km in a 13-day repeating groundtrack. Units are 
signals after PP (right) for recovering hydrology and ice mass variations in cm of equivalent water height 
(top row) and ocean bottom pressure signals (bottom row) given a polar 



Table 3 Signal and error associated with Fig. 4 



Signal (cm) 

Error (cm) 
NoPP 

PP 

Hydrology 

4.67 

13.97 

2.47 

Ice 

2.12 

7.54 

2.06 

Ocean 

2.46 

11.04 

1.71 

Average 

S = 3.08 

E = 10.85 

E = 2.08 


Units are expressed in cm of EWH 


middle and right sets of figures. Furthermore, one can calcu- 
late the power in the truth signals, 5, for hydrology (5(H)), ice 
mass variations (5(1)), and ocean bottom pressure (5(0)) in 
the same manner that the error is calculated in Eq. 1 1 . Table 3 
illustrates the signal and error associated with Fig. 4 in cm 
of EWH for the recovered signals with no post-processing 
(PP) and the recovered signals after PP. In this case, it is seen 
that the error exceeds that of the signal with no PP applied, 
but after destriping and smoothing the solutions the level of 
error is reduced to be lower than the signal. 

Note that throughout the paper, E is obtained by taking 
the solutions to degree 100 and truncating them at degree 60 
to make the spatial maps. The reason for doing this is that it 
was found that if the results are truncated at degree 60, des- 


triping at lower latitudes is generally not required, resulting 
in a gravity solution with better spatial resolution. 


4 Results 

This section shows the most important results from the 
numerical simulations in an effort to optimize k\,ki, and ij. 
Additionally, Sect. 4.5 discusses groundtrack patterns 
obtained by tuning Af2i2 and A m 1 2 • This section begins 
with a discussion on the impact of performing simulations to 
degree 60 versus degree 100. 

4. 1 Degree of estimation 

It was found that the calculated error, E, varies substantially 
for certain cases depending on if the simulations are run to 
degree and order 60 or 100. To explain this, an example case 
is shown. When simulations are carried out to degree and 
order 60, one of the better performing mission architectures 
consists of a polar pair of satellites in an 8-day RP at an 
altitude of 291 km coupled with a lower inclined pair (72°) 
in a 13 -day RP at an altitude of 290 km. The error, E, from 
Eq. 11, is calculated to be 4.64 cm EWH for this case. Con- 
versely, when the simulations are extended to degree and 
order 100, but truncated at degree 60 for a fair comparison, 
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Fig. 5 Logarithm of the error 
in the coefficients for a 
simulation carried out to degree 
60 (left) and a simulation carried 
out to degree 100 (right) 




Order of Clm coefficient 

Fig. 6 Correlations with the C(5 1 ,5 1 ) coefficient for an 8-day RP polar 
orbit 


the error is calculated to be 10.38 cm EWH. This contradicts 
the expected results that the error obtained from these two 
cases should be more or less commensurate with each other. 
To explain this result, Fig. 5 presents the logarithm of the 
actual error in the coefficients as a function of degree and 
order. 

Figure 5 shows that when the solution is extended to 
degree 100, large bands of error show up at two and three 
times the resonant order that do not exist in the degree 60 
solutions. To study why this is the case, the covariance matrix 
can be analyzed. Figure 6 illustrates how the C(51,51) coef- 
ficient is correlated with other coefficients and was obtained 
by examining the covariance matrix of a gravity solution 
involving only the 8-day polar RP orbit over the 13 days 
of the mission simulation. Note that the correlations with 
all coefficients are not shown in the figure, as the correla- 
tions outside of the window shown are effectively zero, as 
expected. 

Given a spherical harmonic coefficient of a certain order, 
coefficients of the same order produce orbital element per- 
turbations of identical frequency (Kaula 1966), shown by 


Table 4 Dominant perturbations for m = 51 and m = 76 for a polar 
pair of satellites in an 8-day RP at 29 1 km 



Period (h) 

A a (cm) 

A® + M (°) 

771 — 5 1 

71 = 51,53,55,... 

7.09 

3.87 

0.86E-5 

p = 24, 25, 26, . . . 
q = 0 
771 = 76 

71 = 77,79,81,... 

7.09 

0.78 

0.67E-6 

p = 36, 37, 38, . . . 
q = 0 





the correlations at order 5 1 in Fig. 6. In this peculiar case, 
however, the period of a near resonant perturbation at order 
m =51 is identical to the period of a near-resonant pertur- 
bation at m = 76, which also manifests itself in Fig. 6 via 
the correlations at m = 76. Table 4 displays the magnitude 
of the largest perturbation in semimajor axis and the along- 
track direction (to + M) for both m = 51 and m = 76 along 
with the period of the perturbation and was calculated via 
Rosborough and Tapley (1987). 

The p and q variables in Table 4 are taken from Kau- 
la’s standard solution for the gravitational potential in terms 
of Keplerian orbital elements (Kaula 1966) while n is the 
degree. The similarity of the perturbation frequencies leads 
to the filter being unable to separate them, and these bands 
of coefficients become poorly determined, as is reflected in 
Fig. 5. This error then manifests itself in the spatial plot of 
the recovered signals, leading to a large value of E. Note that 
similar unexpected correlations exist for orders other than 
the one shown here for this particular case. The above results 
demonstrates the importance of performing simulation stud- 
ies to high degree and order. 

4.2 Selecting an inclination 

It is expected that the selected value for h will have a large 
influence on the mission performance. Figure 7 shows the for- 
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Fig. 7 Covariance analysis for a 17-day RP polar pair coupled with a 17-day RP pair of satellites at various inclinations. The logarithm of the 
standard deviation of the coefficients is shown 


mal errors for a 1 7 -day RP polar pair of satellites coupled with 
a 17-day RP pair of satellites at various inclinations. The log- 
arithm of the formal error of each spherical harmonic coef- 
ficient is shown. Note that the difference in altitude between 
the satellite pairs at different inclinations will have slight 
influences on the formal errors. However, if one examines 
Fig. 3 which shows the altitude for the various inclinations 
of a satellite in a 17-day RP, it is seen that the difference 
in altitude between each of these cases is 15 km at a maxi- 
mum. Thus, this effect should be minimal on the covariance 
analysis. 

Figure 7 shows that for 55° < h < 65°, higher degree 
and order tesseral harmonics are perhaps the best determined. 
Sectorial and near-sectorial coefficients tend to have larger 
errors as h increases. Studying the covariance matrices alone 
might lead one to conclude that an inclination of approxi- 
mately 60° is near-optimal in the sense that the overall errors 
of the coefficients is lowest and the covariance matrix is fairly 
isotropic (no order dependence). However, geographically 
speaking, if the second pair of satellites flies at an inclination 
of 60°, it is seen that it provides no coverage over Green- 
land and does not cover a substantial amount of landmass in 


the northern hemisphere, including Alaska, northern Canada, 
northern Russia, and the Scandinavian countries. Should the 
second pair of satellites provide coverage over these regions, 
it could improve the determination of mass variations in these 
areas even though this is not reflected in the covariance anal- 
ysis. 

Therefore, it is useful to compare the error metric, E, as 
discussed in Sect. 3.2, between different inclinations. Fig- 
ure 8 shows E for a pair of satellites at different inclinations 
in a 17-day RP coupled with a polar pair of satellites in dif- 
ferent repeat periods (15, 16, and 17 days). 

Figure 8 shows that irrespective of whether the lower 
inclined pair of satellites is coupled with a 15-, 16-, or 17- 
day RP polar pair of satellites, the general trend in the error 
as a function of inclination is the same. Typically, the error 
reaches a minimum between 70° and 75°. This is not sur- 
prising since a pair of satellites at this inclination gets fairly 
good coverage over the Earth while still maintaining a sig- 
nificant East- West component in the observable. The results 
from two polar pairs of satellites were not placed on this fig- 
ure as the errors were so high they would distort the scale. 
This attests to the strength of the East-West information in 
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Fig. 8 Error as a function of inclination for a lower inclined pair in a 
17-day RP coupled with a polar pair in 15, 16, and 17-day repeat periods 



Fig. 9 Error as a function of the repeat period of a polar pair of satel- 
lites being coupled with lower inclined pair in a 17-day RP 

the observable. While Fig. 8 shows an extremely small sub- 
set of the results that have been analyzed, the general trend 
of having the best performance for values of n between 70° 
and 75° was consistent across all cases examined. 

4.3 Coupling of repeat periods 

Bender et al. (2008) suggested that having a polar orbiting 
pair in a shorter RP than a lower inclined pair would pro- 
vide more homogeneity in the groundtrack coverage over the 
Earth, and thus, result in better solutions. The polar pair of 
satellites with a shorter RP could reduce the level of temporal 
aliasing errors in the polar regions as well. 

Figure 9 illustrates the error as a function of different 
repeat periods for the polar pair of satellites coupled with 
a lower inclined pair of satellites (70°, 71°, 75°, 80°) in a 


17-day RP. Note that the range of inclinations examined in 
this figure is consistent with the range of inclinations that 
minimized the error as seen in Fig. 8. 

Figure 9 shows that generally the error decreases when the 
lower inclined pair of satellites is coupled with a polar pair of 
satellites in longer repeat periods. Note that this trend exists 
outside of the differences in altitude between the polar pairs 
of satellites. In fact, plotting the error as a function of alti- 
tude of the polar pair rather than the repeat period of the polar 
pair does not show significant trends. It can then be assumed 
that the error is a stronger function of the repeat period of 
the polar pair of satellites than the altitude. Thus, one can 
assume that coupling two satellite pairs with the same repeat 
period will provide near-optimal results. One could argue 
that a global minimum may not be achieved by setting the 
repeat periods of the two pairs equal to each other, based on 
the fact that in Fig. 9 it appears that slightly smaller errors 
exist when the polar pair is either in a 13-day RP or a 16-day 
RP versus the 17-day RP that we have recommended. While 
this is true, the differences in performance between the cases 
is extremely small. In an effort to reduce the search space 
for this type of mission, we feel that invoking a k\ = ki 
constraint leads to near-optimal results while reducing the 
amount of computation time necessary to study all possible 
combinations of k\ and kj. Additional simulation results for 
repeat periods other than 17 days validate this statement (not 
shown). 

4.4 Selecting a repeat period 

After selecting a range of near-optimal inclinations for the 
second pair of satellites as well as enforcing the constraint 
that k\ = ki = L, the search space for an optimal value 
of L is substantially reduced. Table 5 shows the ten cases 
for which results will be displayed. Each of these cases has 
a lower inclined pair with an inclination around 70°- 75°, 
selected to provide the closest altitude to 290 km. Note that 
the case number corresponds to the RP of the satellite pairs 
for convenience. It is evident that there are several repeat 
periods that are not shown. If a particular repeat period is not 
shown, for example, 12 days, this is because the altitude of 
one of the satellite pairs was too high for the results to be 
competitive with those listed. Generally speaking, due to the 
constraint that k/ 1 must be irreducible, repeat periods that 
are prime numbers have a larger range of altitudes to choose 
from. 

Figure 10 shows the error for the cases listed in Table 5 
corresponding to the polar pair of satellites with an altitude 
of Type I. The blue bars are the solutions obtained using 
the processing methodology outlined in Sect. 3.1 (referred 
to as ‘regular’ processing). Comparing these solutions, it is 
seen that having a repeat period in the range of 1 1 to 14 days 
provides the lowest error, with a global minimum provided by 
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Table 5 Mission architectures 
examined to optimize the 
selection of a repeat period 


Case 

Lower inclined 

pair 


Polar pair 




Repeat period 
(days) 

Inclination 

n 

Altitude 

(km) 

Repeat period 
(days) 

Altitude 
( km) Type I 

Altitude 
( km) Type II 

9 

9 

74 

291 

9 

318 

318 

11 

11 

70 

300 

11 

306 

332 

13 

13 

72 

290 

13 

299 

320 

14 

14 

75 

290 

14 

316 

316 

15 

15 

70 

298 

15 

293 

331 

17 

17 

71 

290 

17 

305 

322 

19 

19 

76 

291 

19 

300 

315 

21 

21 

71 

291 

21 

309 

322 

22 

22 

73 

291 

22 

294 

319 

23 

23 

75 

291 

23 

292 

317 



9 11 13 14 15 17 19 21 22 23 

Case Number or Repeat Period 


Fig. 10 Error for the cases listed in Table 5 comparing regular pro- 
cessing and estimating daily 18x18 gravity fields 

L = 13 days. This range of values for L strikes an optimal 
balance between having enough data to form a good solu- 
tion, but a short enough time frame where the accumulation 
of temporal aliasing errors is mitigated. 

The red bars in Fig. 10 are obtained by invoking an alter- 
nate processing methodology where daily 18 x 18 gravity 
fields are estimated in an effort to reduce the level of tem- 
poral aliasing errors. In this process, we simultaneously esti- 
mate degrees 2-18 each day along with degrees 19-100 over 
multiple days (the length of the simulation) during the inver- 
sion process. A multi-day estimate for the lower degree coef- 
ficients is then gained by simply averaging through the daily 
estimates. This process has been shown to reduce temporal 
aliasing errors from high-frequency mass variations at large 
spatial scales and generally results in gravity solutions with 
improved spatial resolution. It is described fully in Wiese 
et al. (201 1). It is seen that estimating the daily gravity fields 


reduces the error substantially for all cases considered. It is 
also interesting to note that with this processing methodology 
invoked, the longer repeat periods provide the lowest errors. 
This makes sense, as for the case with no temporal aliasing 
errors the total amount of error should decrease as the square 
root of the number of observations. However, the reduction 
in errors that a 23-day RP provides over the 1 3-day RP case is 
small when one considers that 10 days of temporal resolution 
are sacrificed. 

4.5 Groundtrack patterns 

The three primary variables in Eq. 8 have been optimized. 
The last two variables, Af2i2, and Azz 12 are not expected to 
have as large of an influence on the solution, as they only 
change the space-time sampling characteristics of the orbit. 
One thing that can be examined, however, is if certain ground- 
track patterns between the two satellite pairs can be devel- 
oped which will lower the errors. 

In this work, it was noticed that for the case of having two 
polar pairs of satellites in the same RP and at the same alti- 
tude, the best solutions are obtained when Af 2 i 2 is set such 
that 

A^12 = S + e, (12) 


where 
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In Eq. 12, 8 shifts the relative node between the two satellite 
pairs exactly 1 80° from each other plus the distance it takes 
for the Earth to rotate during one-half of a satellite revolu- 
tion. The € term is added as a correction factor such that the 
groundtracks of the second pair of satellites will fill in the 
gaps at the equator from the groundtracks of the first pair of 
satellites, resulting in more dense coverage. This architecture 
guarantees that the mutual crossing location of both satellite 
pairs will be at a constant low latitude (~ 7°), rather than at 
the equator. This configuration appears to have substantial 
benefits in the case of two polar pairs of satellites, reduc- 
ing the errors at the resonant orders considerably. The same 
magnitude of improvement is not provided when applying it 
to the case of a polar pair coupled with a lower inclined pair, 
however. The reason for this is twofold: (1) the periods of 
the two satellite pairs are different, and (2) the inclinations 
of the two pairs are different, meaning that the drift rate of 
the node (£2), given in Eq. 15 (neglecting higher order terms), 
is different between the two cases: 

— 3~/jIJ~> R~ cos i 

£2 = — v ~ e 7 . (15) 

2(1 -e 2 ) 2 a 2 

These two differences mean that there are no consistent 
crossings between the two satellite pairs in either space or 
time. However, there are still minor improvements seen when 
invoking Eq. 12 to the cases in Table 5. Figure 1 1 shows the 
reduction in the level of error that this shift provides. Since 
A £2 1 2 will not be constant for the duration of the mission due 
to £2 for the lower inclined pair, the reduction in the level of 
errors seen in Fig. 1 1 represents the natural variability in the 
quality of the solutions due to the precession of A£2n. 

One can now begin to think of developing a spatial ground- 
track pattern for the case of having a polar pair coupled with 
a lower inclined pair that is consistent, as is the case when 



there are two polar pairs. Note that we are only interested 
in developing a spatial pattern (A £2 12 ), and not a temporal 
one (Am 12 ), since the difference in periods between the two 
satellites pairs causes a secular drift rate in this parameter 
(Am 12 ) which cannot be controlled. One can, however, raise 
the altitude of the polar pair of satellites such that its period 
increases enough to compensate for the nodal drift rate of 
the lower inclined pair, thus ensuring consistent crossings at 
the equator in the spatial domain by solving the following 
equation: 

£2ii7ii = <y e (7p P - 71;) (16) 

with 

T = M + w. (17) 

In Eq. 16 and 17, li stands for “lower inclined” and pp 
stands for “polar pair”. The period of the satellites is given 
by T. The constraint given in Eq. 16 guarantees that both 
pairs of satellites complete the same number of orbital rev- 
olutions in the same number of nodal days (/1 = I 2 ) and 
ensures that the groundtracks of the satellites will cross each 
other at constant lines of latitude. Unlike the case of two 
polar pairs, however, the crossings will not have consistency 
in the time domain due to the discrepancy in periods between 
the two satellite pairs. Figure 12 shows the groundtrack of 
two pairs of satellites over South America. The groundtrack 
displayed in blue is from the polar pair of satellites while the 
groundtrack displayed in red is from the lower inclined pair 
of satellites. It can be seen how the two pairs always cross at 
the same latitude. 

The modified architectures necessary to obtain the com- 
plementary groundtrack patterns described above are given 
by the cases involving a polar pair of satellites with an alti- 
tude of Type II in Table 5. Note that the altitude of the polar 
pair for each case, with the exception of Case 9 and Case 14, 
has been raised by approximately 20-30 km with respect to 



Fig. 11 Error for the cases listed in Table 5 comparing the effect of Fig. 12 Complementary groundtrack pattern shown over South Amer- 
shifting the longitude of ascending node ica, arrived at by invoking Eq. 16 


■£) Springer 


94 


D. N. Wiese et al. 


LU 

o 

E 

o 

_c 

m 

o 

LU 


| Regular Cases (Type I) 

| Complementary Groundtracks (Type II) 



9 11 13 14 15 17 19 21 22 23 
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Fig. 13 A comparison of the error between the Type I and Type II 
architectures in Table 5 


Type I. Cases 9 and 14 are the same between Type I and Type 
II by the fact that the polar pair selected was already at the 
appropriate altitude necessary for a complementary ground- 
track pattern. Figure 13 compares the error from the Type I 
architectures in Table 5, with the lowest altitude polar pair 
possible, with those of Type II in Table 5, with a complemen- 
tary groundtrack pattern but a slightly higher altitude for the 
polar pair of satellites. 

Figure 1 3 shows a minor degradation in performance for 
five of the cases considered and a minor improvement in 
performance for three of the cases considered. Cases 9 and 
14 have the same performance since they involve the same 
orbits. It is difficult to draw conclusions from these results. 
Possible benefits from flying the satellites with a complemen- 
tary groundtrack pattern include the fact that the polar pair 
of satellites is at a higher altitude which means increased 
longevity due to lower atmospheric drag forces. Also, the 
crossings at lines of constant latitude could prove beneficial 
in future applications that are not yet realized, e.g., using 
the crossing points as constraint points for determining the 
geopotential at particular locations. For a gravity mapping 
mission, analagous to an altimeter mapping mission such as 
TOPEX or Envisat, there is an argument for having a geome- 
try between the two satellite pairs in terms of the groundtracks 
that permits a consistent synoptic mapping of the time-vari- 
able gravity variations. 

4.6 Expected performance 

While performing an in-depth analysis of the expected scien- 
tific benefits that an optimized architecture consisting of two 
satellite pairs would provide over an architecture consisting 
of a single pair of satellites is out of the scope of this paper, a 
global comparison can be made using the error metrics dis- 


cussed in this paper. For this comparison, we select Case 13 
of Type II in Table 5 consisting of a polar pair of satellites 
coupled with a lower inclined pair of satellites at 72°, both in 
13-day RP orbits and possessing a complementary ground- 
track pattern as discussed in Sect. 4.5. We compare this case 
with one polar pair of satellites in a 13 -day RP at 299 km 
altitude, the performance of which is shown in Sect. 3.2. 
Furthermore, for the case of two satellite pairs, we co-esti- 
mate daily 18 x 18 gravity fields to further reduce temporal 
aliasing errors. As was discussed in Wiese et al. (201 1), one 
benefit of having a polar pair coupled with a lower inclined 
pair of satellites is that temporal aliasing errors can be further 
reduced by co-estimating daily 18x18 gravity fields. Con- 
versely, for the case of a single pair of satellites, it has been 
shown that while co-estimating 10 x 10 gravity fields every 
2 days reduces temporal aliasing errors, the benefits provided 
by making the 2-day estimates are virtually abolished after 
destriping and smoothing the gravity fields. Figure 14 shows 
the truth signals (left), recovered signals (middle), and recov- 
ered signals after PP (right) for recovering hydrology and ice 
mass variations (top) and ocean bottom pressure signals (bot- 
tom) using two satellite pairs. The plots are represented out 
to degree 60 and are expressed in cm of EWH. 

The PP techniques applied to these solutions are different 
than those applied to the single satellite pair solutions. As 
shown in Fig. 14, the addition of the lower inclined pair sub- 
stantially reduces the level of striping in the solutions. There 
are certain bands of coefficients which remain correlated, 
however. The large errors in the high-latitude regions seen in 
Fig. 14 are a direct consequence of correlations in coefficients 
of a fixed order and same parity of degree in the range n > 40 
and 3 < m < 14. Note that these errors occur predomi- 
nantly in geographical areas with latitudes higher than 72° , 
as no East- West information is present here. We found that 
by applying the destriping algorithm presented in Swenson 
and Wahr (2006) to only the selected range of spherical har- 
monic coefficients which have the correlated errors, we are 
able to remove correlated errors at high-latitude regions while 
leaving signals at low latitudes relatively untouched. This 
modified destriping algorithm has been applied to the recov- 
ered hydrology and ice signals in Fig. 14 to represent the 
PP solutions. The PP ocean bottom pressure signals have 
been destriped with this modified algorithm in addition to 
being smoothed with a 200-km averaging radius. Since the 
ocean bottom pressure signals are smaller in magnitude than 
hydrology and ice signals, the signal to error ratio is substan- 
tially smaller and errors manifesting as longitudinal stripes 
are more apparent and hence the need to smooth the solutions. 

Table 6 compares the error from one pair and two pairs 
for the different data reduction methodologies (whether we 
estimate low resolution gravity fields at a high frequency or 
not) along with the different PP techniques. It is seen that 
when no daily (or 2-day) estimates are made, and no PP is 
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Fig. 14 Truth signals ( left), recovered signals (middle), and recovered of Type II in Table 5 after daily 18x18 gravity fields are co-estimated. 

signals after PP (right) for recovering hydrology and ice mass variations Units are in cm of equivalent water height 

(top row) and ocean bottom pressure signals (bottom row) for Case 13 


Table 6 Signal and error 
comparing one pair and two 


Hydrology 

Ice 

Ocean 

Average 

pairs with various data reduction 
methodologies and 

Signal 

4.67 

2.12 

2.46 

S 

= 3.08 

post-processing techniques 

Two pairs error (cm) 







NoPP 

2.55 

5.15 

3.16 

E 

= 3.62 


PP 

2.44 

3.18 

1.44 

E 

= 2.35 


Est.18 x 18 

1.91 

3.31 

2.41 

E 

= 2.54 


Est.18 x 18, PP 

1.87 

1.97 

1.32 

E 

= 1.72 


One pair error (cm) 







NoPP 

13.97 

7.54 

11.04 

E 

= 10.85 


PP 

2.47 

2.06 

1.71 

E 

= 2.08 


Est.10 x 10 

7.80 

4.50 

6.15 

E 

= 6.15 

Units are expressed in cm of 
EWH 

Est.10 x 10, PP 

2.33 

2.24 

1.50 

E 

= 2.02 


performed, the global error, E is decreased by approximately 
67% with the addition of the second pair of satellites. The 
minimal benefits of making 2-day 10x10 gravity field esti- 
mates after PP the solutions for one pair of satellites are 
also seen. Table 6 illustrates that one pair of satellites actu- 
ally has lower errors than two pairs of satellites if only PP 
is performed. The reason for this is because the PP tech- 
niques applied between the two cases are different; if the 
two-pair solutions are destriped and smoothed with a 300- 


km averaging radius, the error metric is actually calculated 
to be E = 1 .99 cm. The solutions with two satellite pairs 
are post-processed (PP) differently in order to preserve the 
information at small spatial scales that would otherwise be 
smoothed over. The best solutions for both cases are obtained 
when daily (or 2-day) estimates are made and the solutions 
have been PP. For this case, it is seen that two pairs of sat- 
ellites offer a more modest reduction in the errors, reducing 
the global error by approximately 15%. Having only a 15% 
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reduction in error can be attributed to the performance metric 
used, which does a poor job of characterizing the distortions 
in the recovered signals from PP. For example, by destriping 
and smoothing (300 km radius) the truth hydrology, ice, and 
ocean signals, we calculate that this PP results in an error of 
E = 1.69 cm. Alternately, applying the modified PP tech- 
niques for two satellite pairs to the truth signals results in an 
error of £ = 0.53 cm, allowing for a more accurate represen- 
tation of the truth signal in the recovery. A visual comparison 
of Fig. 14 and Fig. 4 reveals that the solution from two pairs 
of satellites retains significantly more geophysical signals at 
small spatial scales than the post-processed one-pair solu- 
tions. Note that visually, after PP, the one-pair solutions are 
nearly identical whether we estimate 2-day 10 x 10 grav- 
ity fields or not, with both cases smoothing through signals 
at small spatial scales. As such, it is anticipated that on a 
regional level, one will see much greater benefits with two 
pairs of satellites over one pair of satellites. Future work will 
include an in-depth analysis to quantify the performance of 
adding a second pair of satellites at regional levels for the 
various scientific areas of research. 


5 Conclusions 

Anticipating that future missions dedicated to recovering 
time variable gravity will use laser interferometry for inter- 
satellite ranging, and drag-free technology for altitude con- 
trol, previous studies have shown that temporal aliasing errors 
will dominate the error budget of the mission. One plausible 
way to mitigate these errors is to add an additional pair of 
satellites, increasing the sampling frequency of the mission, 
ultimately leading to a product with greater spatial resolu- 
tion. Additionally, if the second pair of satellites is at a lower 
inclination, the East- West sensitivity of the observable is 
improved, decreasing the longitudinal striping in the solu- 
tions. The goal of this paper was to optimize the orbits of 
two satellites pairs to provide increased spatial resolution in 
determining hydrology, ice mass variations, and ocean bot- 
tom pressure signals globally. 

While the search space for such a problem is, by nature, 
infinite, numerical simulations to degree and order 100 were 
implemented in an effort to reduce it. A search space origi- 
nally consisting of 15 variables was reduced to two variables 
with primary impact on mission performance: the inclina- 
tion of one of the satellite pairs (the other pair is assumed 
to be polar), and the repeat periods of both pairs of satellites 
(shown to be near-optimal when they are equal to each other). 
In this study we considered only circular orbits in repeating 
groundtracks, a minimum allowable altitude of 290 km based 
on a projected 10-year mission lifetime, and assumed a 100- 
km inter-satellite separation distance between each pair of 
satellites. It was found that an optimal value for the inclina- 


tion of the second pair of satellites is between 70° and 75°, 
while an appropriate range for the repeat periods of both 
satellite pairs is between 1 1 and 14 days. The absolute low- 
est errors are given when both satellite pairs are in a 13 -day 
repeat period, one being polar at an altitude of 299 km, and 
the other inclined at 72° at an altitude of 290 km. It should 
be noted that the results of this study are influenced by the 
targeted altitude for the mission as well as the scientific goals 
of the mission. 

The notion of optimizing the relative change in node and 
the argument of latitude between the two pairs was discussed 
in relation to creating complementary groundtrack patterns. 
It was shown that by raising the altitude of the polar pair, 
the nodal drift rate of the lower inclined pair can be compen- 
sated for such that a groundtrack pattern with crossings at 
constant lines of latitude is created. While numerical simula- 
tion results imposing this constraint were not conclusive as to 
whether this definitively results in improved mission perfor- 
mance, there is an argument for having a geometry that per- 
mits consistent global mapping of the gravity field. Finally, 
the importance of extending simulations to high degree and 
order was shown. 

Results showed that with an optimized architecture con- 
sisting of two satellite pairs, global errors are reduced by 
67% (no post-processing) with the addition of the second 
pair of satellites over one pair of satellites. After each set of 
solutions has been destriped and smoothed with appropriate 
methods, and daily (or 2-day) gravity fields are estimated, it 
was seen that two pairs of satellites offered a more modest 
reduction in error of 15% over one pair of satellites. Visually, 
however, it was seen that the two-pair solution retains sig- 
nificant geophysical signals at small spatial scales which are 
smoothed and damped in the one-pair solutions. Future work 
will involve an in-depth examination of the expected scien- 
tific benefits of an optimized two-pair architecture, extending 
the analysis to local regions as well as longer time spans. 
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